#include "mask.h"
#include "pattern.h"

double dc(double l, double w, double g, double d, double r, double o)
{ //draw a dc coupler with I/O s-bend
  //l: coupling lenth
  //w: waveguide width
  //g: coupler gap
  //d: s-bend offset
  //r: s-bend bent radius
  //o: s-bend offset
  
  double dy;
  double dx;
  push_cp();
  dy = (d*2+w+g)*0.5;
  // upper line
  offset(-dy);
  dx = sbend(d, r, o, w);
  sw(l, w);
  sbend(-d, r, o, w);
  pop_cp();
  
  //lower line
  offset(dy);
  sbend(-d, r, o, w);
  sw(l, w);
  sbend(d, r, o, w);
  
  offset(-dy); //leave the cp at the center of the two outputs
  return dx*2+l;//return the x span
}

double sbend2(double d1, double d2, double r, double o, double w)
{ //draw two s-bend with distance from d1 to d2

  double z;
  push_cp();
  offset(-d1*0.5);
  z = sbend((d1-d2)*0.5, r, o, w);
  pop_cp();
  offset(d1*0.5);
  sbend(-(d1-d2)*0.5, r, o, w);
  offset(-d2*0.5);
  return z;
}

void taper_sbend_v1(double l, double w1, double w2)
{ //taper with the sbend edge.
  const double dx = getlayer_accuracy(currentlayer());
  if (fabs(w1-w2) < 4*dx || fabs(l) < 2*dx )
  {
    taper(l, w1, w2);
    return;
  }
  if (w1>w2){
    skip(l);
    push_cp();
    rotate(M_PI);
    taper_sbend(l, w2, w1);
    pop_cp();
  } else {
    double y1 = (w2-w1)*0.5*0.25;
    double x1 = l*0.25;
    double r = y1 + x1*x1/y1;
    
    y1 = w1;
    double y2;
    
    x1 = 0;
    double x2=dx;
    while (x2<0.5*l)
    {
      y2 = w1+2*(r-sqrt(r*r-x2*x2));
      //if ((y2-y1)>2*dx)
      if ((y2-y1)>4*dx+4*sqrt(dx*(y1+dx)))
      { x2 -= dx;
        y2 = w1+2*(r-sqrt(r*r-x2*x2));
        taper(x2-x1, y1, y2);
        
        y1 = y2;
        x1 = x2;
      }
      x2 +=dx;
    }
    y2 = 0.5*(w1+w2);
    taper(0.5*l-x1, y1, y2);
    
    x1 = 0.5*l;
    y1 = 0.5*(w1+w2);
    x2 = x1 - dx;
    while (x2>0)
    {
      y2 = w2-2*(r-sqrt(r*r-x2*x2));
      //if ((y2-y1)>2*dx)
      if ((y2-y1)>4*dx+4*sqrt(dx*(y1+dx)))
      {
        x2 += dx;
        y2 = w2-2*(r-sqrt(r*r-x2*x2));
        taper(x1-x2, y1, y2);
        
        y1 = y2;
        x1 = x2;
      }
      x2 -=dx;
    }
    y2 = w2;
    taper(x1, y1, y2);
  }  
}

void taper_sbend(double l, double w1, double w2)
{ //taper with the sbend edge. still need to be improved by using low-level function. 
  const double dx = getlayer_accuracy(currentlayer());
  if (fabs(w1-w2) < 4*dx || fabs(l) < 2*dx )
  {
    taper(l, w1, w2);
    return;
  }
  if (w1>w2){
    skip(l);
    push_cp();
    rotate(M_PI);
    taper_sbend(l, w2, w1);
    pop_cp();
  } 
  else 
  {
    double y1 = (w2-w1)*0.5*0.25;
    double x1 = l*0.25;
    double r = y1 + x1*x1/y1;
    double a = 2*atan(y1/x1);    
    double da = 2*acos(r/(r+dx));
    int nstep = (int)(fabs(a)/da)+2;
    da = a/(nstep-1);
    
    double dl, y2;
    y1 = w1;
    for(int i=1; i<nstep; i++)
    { 
      dl = r*(sin(i*da)-sin((i-1)*da));
      y2 = (r-r*cos(i*da))*2+w1;
      taper(dl, y1, y2);
      y1 = y2;
    }
    //fprintf(stderr, "%e, %e, %e, %e", y1, (w1+w2)*0.5, r, a);
    y1 = w1*0.5 + w2*0.5;
    for(int i=nstep-1; i>1; i--)
    { 
      dl = r*(sin((i)*da)-sin((i-1)*da));
      y2 = w2 - (r-r*cos((i-1)*da))*2;
      taper(dl, y1, y2);
      y1 = y2;
    }        
  }  
}

void pie(double r, double a, double a0)
{//draw a pie with initial angle of a0, cp at the circle center
  push_cp();
  rotate(-a0);
  skip(0.5*r);
  rotate(((a>0)?0.5*M_PI:-0.5*M_PI));
  cw(0.5*r, fabs(a), r);
  pop_cp();
}

void star_coupler(int Nin, int Nout, double Ry, double dx, double D, double w, double w2, double l, bool flag_horn)
{//Nin x Nout star coupler with tapered I/O.
  double dain = dx/Ry;
  double daout = D/Ry;  
  double arc_out = D*(Nout-1 + 2);
  double arc_in = dx*(Nin-1 + 2);
  
  push_cp();
  skip(l);
  
  push_cp();
  rotate(-arc_out*0.5/Ry);
  skip(0.5*Ry);
  rotate(0.5*M_PI);
  cw(0.5*Ry, arc_out/Ry, Ry);
  pop_cp();
  
  push_cp();
  rotate(-(D*(Nout-1)/Ry*0.5));
  for(int i=0; i<Nout; i++)
  {
    push_cp();
    skip(Ry-l);
    if (flag_horn)
    {
      taper_parabolic(2*l, sqrt(2*w2*w2-w*w), w);
    }
    else
    {
      taper(2*l, 2*w2-w, w);
    }
    pop_cp();
    rotate(daout);
  }
  pop_cp();
  
  skip(Ry);
  rotate(-M_PI);
  
  push_cp();
  rotate(-(dx*(Nin-1)/(Ry)*0.5));
  for(int i=0; i<Nin; i++)
  {
    push_cp();
    skip(Ry-l);
    if (flag_horn)
    {
      taper_parabolic(2*l, sqrt(2*w2*w2-w*w), w);
    }
    else
    {
      taper(2*l, 2*w2-w, w);
    }
    pop_cp();
    rotate(dain);
  }
  pop_cp(); 
  
  push_cp();
  skip(Ry*0.5);
  rotate(-arc_in*0.5/(0.5*Ry));
  skip(0.25*Ry);
  rotate(0.5*M_PI);
  cw(0.25*Ry, arc_in/(0.5*Ry), Ry*0.5);
  pop_cp();
  
  double l1 = Ry*cos(arc_out/Ry/2);
  double wt1 = 2*Ry*sin(arc_out/Ry/2);
  skip(Ry-l1);
  double l2 = 0.5*Ry*cos(arc_in/(0.5*Ry)/2);
  double wt2 = 2*0.5*Ry*sin(arc_in/(0.5*Ry)/2);
  taper(l1-(Ry*0.5-l2), wt1, wt2);
  
  pop_cp();
  skip(l*2+Ry);
  
}

/*
void AWG(int Nin, int Nout, double lambda0, double s_GHz, double dx, double D, double w, double ns, double Ns, double nw, double Nwg)
{
  double clight = 299792458.;
  double Ry = ns* dx* D* Nout/lambda0;
  double dL = ns*clight/lambda0*dx *D/(Nwg*s_GHz*Ry);
  double m = nw * dL/lambda0
}*/
